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Abstract. We describe a fast method to simulate in a semi-analytical 
way consistent maps of the thermal, kinetic and polarised Sunyaev-Zel'- 
dovich effect, featuring both cluster spatial correlations and large-scale 
velocity flows. 



1. Introduction 

There are a large number of applications for simulated maps of the Sunyaev- 
Zel'dovich (SZ) emission of galaxy clusters. Such maps are being used, for 
instance, to assess the capabilities of instruments dedicated to the SZ observa- 
tions, or for the development of data analysis tools. They can also be used to 
constrain models of structure formation and/or cosmological parameter sets by 
comparison to real data. 

There already exists a large number of tools to simulate SZ maps (Bond et 
al. 1996, Aghanim et al. 1997 , da Silva et al. 2000, Kay et al. 2001, Kneissl et al. 
2001, da Silva et al 2001). Such tools range from simple uncorrelated Poisson 
distributions of idealised clusters to sophisticated hydro-codes simulating the 
complex physical process of structure formation, including non-linear gravity 
and some additional ingredients, such as magnetic fields, dissipation. . . While 
the former may be over simplistic for a range of applications, the latter require 
large computing facilities and extensive computer time, and do not illuminate 
simply the influence of the code inputs on the final results. In the approach 
presented here, we develop a code for simulating SZ maps that is fast enough 
to permit numerous Monte-Carlo simulations in a wide range of cosmological 
scenarios, while being accurate enough to correctly account for all three SZ 
effects (thermal, kinetic and polarised) and their distribution in both real and 
velocity space (e.g., bulk flows) according to the underlying spectrum of density 
fluctuations P(k). 



2. The SZ effect and galaxy clusters 

Observations of galaxy clusters are useful as a probe of cosmology, and obser- 
vations of the SZ effect in particular present various advantages in this light. 
Most notably, as galaxy clusters are the largest bound objects in our Universe, 
their distribution as a function of redshift is extremely sensitive to the history of 
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large-scale structure formation. With its distance independence, the SZ effect 
is ideal for tracing this history out to large redshifts, and thereby for providing 
stringent constraints on models. The SZ effect, although for long below attain- 
able instrumental sensitivity, now represents an established and unbiased tool for 
the detection of distant clusters and for the estimation of their hot gas content. 
In particular, the Planck mission, to be launched by ESA in 2007, will produce 
an all-sky catalog of several tens of thousands of clusters, which will be used 
to constrain scenarios of structure formation and cosmological parameters, for 
example, through cluster count tests. In addition, SZ-selected cluster catalogs 
will be useful to constrain the two-point correlation function of clusters, and 
hence the matter power spectrum P(k) on the largest scales. 

The kinetic and polarized SZ effect, although much more difficult to detect 
than the thermal SZ, also appear as extremely promising tools, allowing the 
determination of both radial and transverse velocities for clusters. Although 
detecting the polarised SZ effect still appears only marginally possible (Audit 
& Simmons 1999, Sazonov h Sunyaev 1999), the measurement of bulk radial 
velocities through the kinetic SZ effect seems within the reach of upcoming 
instruments (Aghanim et al. 1997). Both will become of ever greater importance 
for cosmology as instruments improve. 

2.1. Why a new code for simulating SZ effect maps 

A key issue in interpreting SZ observations is the understanding of the selec- 
tion function for clusters and of the importance of measurement errors. The 
assessment both requires testing methods for the extraction of the relevant SZ 
observables on simulated observations that accurately model the cluster popu- 
lation - gas physics and distribution in both real and velocity space - as well 
as instrumental effects and relevant astrophysical foregrounds. The simulations 
should be performed over a wide range of cosmological scenarios. As an example, 
one expects that in the near future direct constraints of cosmological parameters 
will be obtained from SZ cluster counts and redshift distributions. Likelihood 
contours may be obtained using analytical predictions of cluster counts (or clus- 
ter distributions in the observable parameter space), but the approach must 
be tested using Monte-Carlo simulations with numerous realisations of the ob- 
servations and subsequent data analysis. Clearly, this can not be done using 
full-scale structure formation codes simulating gravitational collapse, as such 
codes require days or weeks of computation on supercomputers. 

Our goal is to simulate maps accurate enough that they reliably feature the 
physical properties that convey cosmological information (e.g., cluster counts, 
cluster correlations in real space, bulk flows), in computation times permitting 
one to test methods statistically on thousands of simulations over a wide range 
of cosmologies. 

2.2. The formation of cluster sized dark-matter haloes 

The formation of large haloes takes place through the gravitational instability 
of dark matter density fluctuations. Original density fluctuations (at the time of 
matter radiation decoupling), still well within the linear regime, grow through 
gravitational instability into the non-linear regime (when the density contrast 
5 = 5p/p becomes of order unity). Collapsed haloes of dark matter appear 
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where gravitational attraction overcomes the expansion and the density contrast 
becomes larger than some limit. The Press-Schechter (PS) formalism (1974) 
predicts the average number of such dark haloes as a function of mass, M, and 
redshift, z, for Gaussian fluctuations: 
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(1) 



where v = S c (z)/a(z, M), p is the present average matter density, dV c the 
covolume element, 5 c (z) ~ 1.68 the critical density contrast necessary for the 
collapse, and a(z, M) = D g (z)a(0, M) is the linear density fluctuation amplitude 
evolving according to the linear growth factor D g (z) (= 1 at z = 0). The 
fluctuation amplitude a(z, M) is related to the power spectrum P(k) by: 
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where W(k) is the Fourier transform of a 3D tophat smoothing window func- 
tion of size R = (3M/47T/?) 1 / 3 . The power spectrum P(k), however, enters 
this formalism only through its integral value (the variance of the fluctuations); 
therefore, the formalism cannot account for the spatial correlation of clusters. 
Simulations based simply on PS result in maps that can only reproduce the pre- 
dicted number of collapsed halos as a function of redshift, and not the cluster 
two-point correlation compatible with the underlying P(k). 

The power spectrum uniquely defines the two-point correlation properties 
of the density perturbation field, S. Kaiser (1984) was the first to suggest that 
clusters trace the density field in a biased fashion, preferentially forming in 
over dense regions: 

^clusters OC 1 + b5 (3) 

where the bias parameter b(M, z) is in general a function of both mass and 
redshift. This behavior is explained by the fact that galaxy clusters form from 
only the rarest, highest density peaks in the density field, and it is borne out 
by large N-body simulations; for example, Sheth et al. (2001) give an analytic 
expression for b(M, z), based on earlier work of Mo & White (1996), that matches 
the clustering seen in simulations. We note in passing that, also from studies 
of clustering in large N-body simulations, the Press-Schechter formula is known 
to somewhat over-predict the number density of small mass dark haloes and to 
underestimate large haloes. These studies have lead to improved expressions for 
halo number density (e.g., Jenkins et al. 2001). 

2.3. Cluster velocities 

Cluster velocities are needed for the simulation of both kinetic SZ (radial part) 
and polarised SZ maps (transverse part). In the simplest simulation approach, 
the velocity of each individual cluster can be drawn at random from a 3D Gaus- 
sian probability distribution with rms proportional to the square root of the 
integrated velocity power spectrum: 



a v (z,M)=a(z)H(z)f(n m ,A) 
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Figure 1. Thermal SZ map of 3 by 4 square degrees for a universe 
with Q m = 0.3, A = 0.7, h = 0.65, tt b = 0.05 (left panel), and with 
£l m = 1, A = 0, h = 0.5, Qb = 0.07 (right panel). Note that because 
of the different baryon fraction, the A universe clusters are brighter 
in SZ (respective color scales have been chosen such that they range 
from y = 0toy = 4x 10~ 5 on the left panel , and from y = to 
y = 1.5 x 10~ 5 on the right one, with higher values saturated. 

with a(z) = l/(l + z) and H(z) = #oV^m(l + z f + (1 - - A)(l + z) 2 + A. 
The function / is defined as / = din D g {z) / dlna. Such an approach results in 
velocities which have the proper magnitude, but which are not correlated from 
one cluster to the next (i.e., no bulk flows). For sensitivity reasons, and because 
of confusion due to the primordial CMB fluctuations (and other foregrounds), 
upcoming experiments will not be able to measure individual cluster velocities 
(Haehnelt and Tegmark 1996, Aghanim et al. 2001). An estimate of the ability 
of instruments like Planck to measure peculiar velocities on large scales can 
be obtained by extrapolating the expected accuracy on a single cluster to a 
larger number of clusters moving together in a bulk flow (Aghanim et al. 2001). 
However, a more accurate study of this approach and its uncertainties requires 
simulations that correctly model large-scale velocity correlations. We detail 
such a method in the next section. 

3. Simulation method 

3.1. The cluster spatial distribution 

Since the power spectrum of density fluctuations, P(k), uniquely defines the 
density field responsible for the formation of halos and their clustering, we start 
with a random realisation of the density contrast in a 3D comoving box, with the 
user placed at one end. The size of the box is matched to the size of the required 
map (angular extent). For a 3 degree-by-3 degree map, typical comoving box 
sizes (extending to a redshift of 5) are 600 x 600 x 6000 Mpc, divided in resolution 
cells of about 20 Mpc. To account for the different redshifts seen by the observer 
while looking through the box, the density field is scaled by the linear growth 
factor, D g (z), over slices of redshift (distance from the user). A random catalog 
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within the box is first constructed by drawing from a Poisson distribution with 
average given by the PS law (or Sheth-Tormen, or Jenkins et al., as appropriate, 
up to the user of the simulation tool) the number of clusters in a grid of mass 
bins for each redshift shell. The average can be slightly modified to take into 
account the fact that the particular redshift shell might be slightly overdense 
or underdense with respect to the ensemble average. For each shell, we then 
distribute the clusters in (9, <j)) with a probability proportional to 1 + b5, where 
S is the density field scaled to the redshift of the shell and b the clustering bias. 
In this way, we construct a cluster distribution compatible with the underlying 
P{k) (within the resolution of the code). 

3.2. Cluster peculiar velocities 

The peculiar velocity of each cluster is found using linear theory applied to the 
density field realisation, 5p/p. With the assumption that the peculiar velocity 
field is irrotational, as expected in most scenarios of structure formation, the 
velocity field can then be obtained from the Newtonian gravitational potential 
through 

3n m (z)H(z) a ( ' 

where the gravitational potential $ is connected to 5 = 5p/p by 

V 2 $ = ^Ga 2 5p (6) 



3.3. Cluster gas properties 

The physical modeling of individual clusters is a complicated issue, in principle 
calling into play a variety of physical processes. As spectacularly illustrated by 
the new X-ray satellites XMM-Newton and Chandra, clusters are not simple 
objects, but rather display complex substructures and inhomogeneities. For 
large maps which will be observed at low angular resolution (as will be the case 
for the Planck mission, for instance), these substructures are lost and thus not 
very important for the final properties of the map. For resolved clusters, on the 
other hand, proper physical modeling of cluster substructure is an important 
issue. 

In our current simulations, we adopt a simple, self-similar cluster model 
based on an isothermal (3 model: 

/ / \ 2\ -W/ 2 

n e (r)=n e (0) 1+ - (7) 



for the density profile, and a gas temperature-mass relation normalized to nu- 
merical simulations (Evrard et al. 1996): 

where Anl(z) is the non-linear density contrast on collapse (a weak function of 
both O and A that equals 178 in a critical universe). The virial radius is given 
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Figure 2. Velocity map of 3 by 4 square degrees. 
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by (e.g., Bartlett 1997) 



_2/ 3 / M \ 1/3 / 178 X 1 / 3 



r„ = 1.69/T 2/d TT — (1 + z)" 1 — — Mpc (9) 

Another possible approach is to replace this simple cluster model by maps of 
the Comptonization parameter y and / or the optical depth r obtained in detailed 
simulations of individual clusters, as done by Kneissl et al. (2001). 



4. Examples of maps obtained with the software 

In Figure 1 we show 3x4 degree 2 SZ maps simulated with our code for two 
different cosmologies, one with (fi m = 0.3, A = 0.7, h = 0.65, f2& = 0.05), 
left panel, and one with (£l m = 1, A = 0, h = 0.5, = 0.07), right panel. 
The visual difference between the two maps is striking, with more high redshift 
clusters in the A-model than in the critical model (Barbosa et al. 1996). 

As another example, we show in Figure 2 a peculiar velocity map where 
radial velocities are coded in color (the map is one of radial velocity times optical 
depth, (3 r r). Transverse velocities for clusters between redshifts 0.8 and 1 are 
indicated by the small arrows starting from the center of each cluster. The 
tangential velocities in this limited redshift slice demonstrate clear correlations 
on large scales. Maps of the polarised emission are straightforwardly obtained 
from the optical depths and transverse velocities, but are not shown here due to 
a lack of space. 



5. Characteristics of the code 



The software can produce simulations for critical, open or closed CDM-like 
models, with or without a cosmological constant. A future development will 
be the inclusion of an arbitrary equation-of-state for the dark energy term 
(Quintessence). Presently, our software can only produce small tangential maps 
of limited fractions of the sky, although we are currently working on full-sky 
extensions. The software is implemented in the IDL programming language, 
and requires only a couple minutes on a modest PC to produce a map of several 
degrees squared. 



6. Conclusion 



We have written a fast code for the simulation of SZ maps (thermal, kinetic 
and polarised) based on a simple model of cluster physics and on Monte Carlo 
realizations of the mass function and the linear density field. The code produces 
a cluster distribution and large-scale peculiar velocities consistent with the un- 
derlying power spectrum, P(k). The software can produce maps for any kind of 
open, closed or critical universe, with or without a cosmological constant. With 
the expectation of data from the new generation of SZ dedicated instruments, 
this will be a useful tool for a wide range of applications, from optimizing obser- 
vations strategies and devising/testing data processing methods, to constraining 
cosmological scenarios with real data. 
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